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Abstract 

We compute the spectral statistics of the sum H of two independent complex Wishart matrices, 
each of which is correlated with a different covariance matrix. Random matrix theory enjoys 
many applications including sums and products of random matrices. Typically ensembles with 
correlations among the matrix elements are much more difficult to solve. Using a combination of 
supersymmetry, superbosonisation and bi-orthogonal functions we are able to determine all spectral 
/c-point density correlation functions of H for arbitrary matrix size N. In the half-degenerate case, 
when one of the covariance matrices is proportional to the identity, the recent results by Kumar for 
the joint eigenvalue distribution of H serve as our starting point. In this case the ensemble has a 
bi-orthogonal structure and we explicitly determine its kernel, providing its exact solution for finite 
N. The kernel follows from computing the expectation value of a single characteristic polynomial. 
In the general non-degenerate case the generating function for the fc-point resolvent is determined 
from a supersymmetric evaluation of the expectation value of k ratios of characteristic polynomials. 
Numerical simulations illustrate our findings for the spectral density at finite N and we also give 
indications how to do the asymptotic large-A analysis. 
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1 Introduction 


One of the main goals and achievements of random matrix theory is the quantification of noise in the 
spectral statistics of some given data or operator. It allows to make analytical predictions for typical 
fluctuations and was called “a new kind of statistical mechanics” in [T] . One of the ensembles of random 
matrices most studied in physics and mathematics was introduced by Wishart [2] in the context of 
mathematical statistics and is frequently used in multivariate statistics [3]. Probably the best known 
quantity of the Wishart ensemble is its global spectral density derived by Marchenko and Pastur [1], 
which is obtained in the limit of infinite matrix size. This and other more refined quantities like 
eigenvalue correlation functions or distributions of individual eigenvalues have found applications in 
physics, e.g. in Quantum Chromodynamics [5] and in other disciplines such as hnance [6], medicine (7j 
climate research [8], telecommunication [9], cf. |10] for a comprehensive list of modern applications. 

In the classical Wishart ensemble a minimum of information is imposed by assuming that all its 
matrix elements are independent real or complex normal random variables. Due to this property this 
ensemble is exactly solvable for finite matrix size, see e.g. m- However, in realistic data system 
specihc correlations among matrix elements are observed. In order to implement these the correlated 
Wishart ensemble has been introduced [3]. Here the eigenvalues and eigenvectors become coupled in 
a non-trivial way. Real matrices pose a considerable challenge as the group integrals involved in the 
calculation are not available. Nevertheless the spectral density [12] has been computed exactly in this 
setting for finite matrices. For complex matrices the exact solvability persists, see |13j . 

Further generalisations that encode more structure of the observed correlations quickly lead to 
more difficult random one- or multi-matrix models. This is why we stick to complex matrices. The 
model introduced by Kumar m that we will study and solve in the present work belongs to this 
class. Let us now explain under which assumptions it arises. Consider an x Nw matrix Wgt 
that contains N different time series in its rows, measured at Ny/ times steps given by its columns. 
When averaging over different realisations as denoted by brackets (...) the most general correlations 
between two time steps t and t' and two time series s and s' read {WstWs’t') = Several 

approximations to this situation have been proposed. First, it was assumed that correlations in 
time steps and among time series factorise, as was considered in economics |15] . climate research |8|, 
sociology |16] and telecommunication |9]. In random matrix theory this problem was analytically 
discussed recently in [EldHI also for the real case and is called doubly-correlated Wishart ensemble. 
Second, a further strong simplihcation is to assume that different times steps are uncorrelated. Here 
analytical derivations for the spectral density were obtained in naiiiiiisi. In this simplihed case 
also cross-correlations were considered |20ll21j . in particular when the correlations among time series 
exhibit a matrix block-form, see also m for related time-lagged correlations. 

In a very recent work by Kumar m the assumed absence of correlations among different time 
steps was softened. Two epochs of Na and Nb time steps with Na + Nb = Ny/ were introduced, 
with different spatial correlations during the two epochs, see also |23] for an earlier work where such 
ensembles occur. This is equivalent to study the statistics of an iV x iV random matrix H which is given 
by the sum of two correlated Wishart matrices, or for more epochs T >2hy H = For two 

epochs, T = 2, Kumar determined the joint density of eigenvalues of H in terms of a hypergeometric 
function of matrix argument. In the case when one of the two correlation matrices is proportional 
to the identity matrix, called half-degenerate case, he showed that this joint density reduces to a bi- 
orthogonal ensemble containing ordinary confluent hypergeometric functions of Kummer type. In this 
case he gave determinantal expressions of size N -\-k containing moments of Kummer’s hypergeometric 
functions for density with k = 1 m and for the fc-point eigenvalue correlation functions [24|. Only 
in the completely degenerate case of two equal covariance matrices H reduces to a single rectangular 
correlated Wishart matrix, where the spectral statistics are well known. 
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Our goal in the present work is to first exploit the bi-orthogonal integrable structure of this 
ensemble at half-degeneracy in the sense of [25], which was only mentioned in M- By explicitly 
determining the kernel of bi-orthogonal functions we provide the solution for all fc-point correlation 
functions in terms of a determinant of size k for arbitrary N. In this form the exact solution is 
amenable to study the asymptotic large-A^ limit. Second, we aim at the solution for the general 
non-degenerate case that was not considered in |141I24| . Applying supersymmetric techniques [261128] 
together with superbosonisation, we provide a closed form for the generating function of the fe-point 
resolvent. 

One particular application, where the solution of this ensemble is relevant for hnite but not nec¬ 
essarily small N, is in the area of wireless telecommunication, in so-called multiple-input-multiple- 
output systems (MIMO) |29It31j . In the setup of T = 2 the sum of two correlated Wishart matrices 
corresponds to the situation of two sets of transmitters and one set of receivers, with an obvious 
generalisation for larger T. The central quantity studied in MIMO systems is the mutual information 
averaged over different realisations called ergodic capacity. This quantity can be expressed as the 
expectation value of In [det [1 -|- H]]. Its generating function is closely related to those of the fe-point 
resolvents that we compute. 

Our work is organised as follows. In section [2] the ensemble we study is introduced and previous 
results from m and [231 are summarised. In section [3] we present our solution for the half-degenerate 
case. This includes a closed form expression for the partition function derived in appendix [A] for the 
kernel in subsection [Q and for the polynomials orthogonal to Kummer’s hypergeometric function 
in subsection 13.21 These polynomials are given by the expectation value of a single characteristic 
polynomial. We determine the latter for the general non-degenerate case. The expectation value of an 
inverse characteristic polynomial follows in subsection 13.31 The spectral density in the half-degenerate 
case is given in detail in subsection 13.41 and illustrated together with Monte Carlo simulations. The 
generating function for the /c-point resolvent in the non-degenerate case is presented in section U 
with details described in appendices [C] and |Dj As an example we give an explicit representation for 
the spectral density and sketch its large-A^ asymptotic analysis. At the end of section |3] further 
generalisations are briefly discussed. In section [5] we summarise our results and present an outlook. 

2 Formulation of the Problem 

In this section we dehne the random matrix ensemble that we consider and summarise the results 
derived by Kumar [141124] . which serve as our starting point. Our notation follows closely the one 
employed in [14] . 

We begin with two independent copies of complex correlated Wishart ensembles composed of 
Gaussian random matrices A and B of dimensions N x Na and N xNb, respectively. Their dimensions 
are restricted to Na > N and Nb > N, and they are distributed as 

Va{A) = , 

Vb{B) = det[SB]“^''e"^^(^s'^^^) . (2.1) 

The probability distributions (12.ip represent two independent complex correlated Wishart ensembles 
AA^ and BB\ respectively. Furthermore we require throughout our discussion that the fixed corre¬ 
lation matrices Sa ^ are positive dehnite, which should be by dehnition the case. The partition 
function is normalised to unity and it completely factorises, 

2n = J [dA] j [dB] Va {A) Vb {B) = l. (2.2) 
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Integration measures over random matrices denoted by [dX] are the Lebesgue measure meaning the 
product over the differentials of all independent real and imaginary parts of the elements of X. 

We are interested in the spectral statistics of the sum of the two independent Wishart matrices, 

H = AA^ + BB^ = WW\ (2.3) 

which is X Hermitian and positive definite. Here we have introduced the larger rectangular matrix 
W = {A,B) of size N x {Na + Nb)- It consists of the elements Wij = Aij and = ^i,k for 

l<z<A^, l<j< Na and 1 < k < Nb- This definition will be helpful later on due to the duality of 
the matrices WW^ and W^W, meaning that both matrices share the same non-zero eigenvalues. 

In [14] two equivalent representations for the distribution of matrix elements of ff were shown to 
hold in the general case, valid for arbitrary non-degenerate correlation matrices T^a and 

Vh{H) = Ch det[H]”*e-Tr (^Nb-,Na + Nb; H) 

= Ch det[H]-e-'^ {Na; Na + Nb; {^b^ - H) . (2.4) 

They are given in terms of the confluent hypergeometric function iXi of matrix argument, cf. [32] for 


a definition, together with 



m 

= Na + Nb-N, 

(2.5) 

^-1 

N 

= detpA]^^ detpse]^"" JJ F(m + j). 

_ 1 

(2.6) 


i=i 


For later use let us also define expectation values in the general non-degenerate case, 

= / [dH]0{H)VH{H ). (2.7) 

Here, 0{H) is a function of the matrix H, e.g. the characteristic polynomial det [xItv ~ H] that plays 
an important role later. 

In the half-degenerate case where one of the covariance matrices is proportional to the identity, 
i.e. Tja = and = diag((TBi,... ,aBN), it was shown in [T3| that the hypergeometric function 

in eq. ()2.4p reduces to a determinant of Hummer’s confluent hypergeometric functions iFi, using an 
identity from [33|. The joint probability density of eigenvalues Xj, j = 1,... ,N, of H can then be 
written as 

N 

Pn{Xi,...,Xn) = AnUXi}) 

i=i 

xdet [iFi {m + 1 - Na; m + 1; \i<k,i<N] - (2.8) 

Furthermore, in m an expression for iFi at the above parameter values was given as a sum over 

incomplete Gamma-functions. In eq. (12.8p we also have introduced the Vandermonde determinant 

AAr({Ai}) = (Aj — Aj) = det[A:^ (2-9) 

Following m, the normalisation constant in eq. (12.8p is given in terms of a determinant of Gauss’ 

hypergeometric function 2 F 1 , 

‘'A ’ + ( 2 . 10 ) 

k=l 

X det[2Ti {m + 1- NA,m + l+j;m + l; - cfb])(ta) \i<i,j<N] ■ 
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In the completely degenerate limit asj —>• cta for all j, the differences — <^Bi) and thus both 
determinants in the normalised joint probability density, eq. (|2.8p with eq. (12.101) . vanish. L’Hopital’s 
rule eventually reduces the limiting expression to the joint probability density of a single uncorrelated 
Wishart ensemble, which is also called Laguerre or chiral Gaussian unitary ensemble. 


N 

lim PAr(Ai,..., Aat) ~ IT AAr({Aj})^ . (2-11) 

j=i 

This classical ensemble can be solved in terms of generalised Laguerre polynomials as orthogonal 
polynomials, thus its name. 

More generally it can be easily seen that in the limiting case of equal correlation matrices, = 
Ss = S, not necessarily proportional to the identity matrix, the joint probability density of H reduces 
to 

Vh{H) ~ det[i/]™e"'^ , (2.12) 

corresponding to a single rectangular correlated Wishart ensemble. 

We underline that the completely different limit, where —)• ctb / for all j, does not yield 

the Laguerre ensemble as immediately follows from the representation (|A.14I) of Kummer’s confluent 
hypergeometric function. Then, the one-point weights resulting from eq. (j2.ip consist of two different 
exponentials, namely exp[—cr^^^Aj] and exp[—instead of one. 

As it was mentioned in m the joint probability density (12.8h in the half-degenerate case represents 
a bi-orthogonal ensemble in the sense of Borodin [25]. Consequently it satisfies a determinantal point 
process. Following Mehta [34| the /c-point density correlation functions are defined by 

Nl A 

Rk (Ai,..., Afc) = _ / dXjPN{Xi,---,XN) (2-13) 

1 j=k+i'^^ 


and can be expressed in terms of a single kernel |25|, once the bi-orthogonal functions are determined. 
For the proof of existence and equivalence of general systems of bi-orthogonal functions see [35| . We 
will pursue this approach in section [3] below. 

Following a different approach, in [TTj in the half-degenerate case the spectral density i?i(Ai) of 
H was expressed as a determinant of an (N -|- 1) x (N + 1) matrix by using the generalisation of the 
Andreief integral formula, see [36] or a summarised version in our appendix [Bj With the help of the 
same formula this was then generalised in [23] to the fc-point density correlation function expressing 
it as a determinant of an (N + k) x (N + k) matrix: 


(Ai,..., Xk) = det 


/o°°dAA^-VTA)| 


(Ai)| 







Ofcx/c 


(2.14) 


where 


(pj{X) = A™e iFi (m + 1- NA]rn + l] (cj^^ - cr^j)A^ 


(2.15) 


with the normalisation constant given by ()2.inp . Whilst these representations are valid expressions, 
that may be useful for small matrix size N, they are clearly not suitable to take the large-Al limit 
or to study the issue of universality. In particular they do not exploit the integrable structure of a 
bi-orthogonal ensemble [25] , expressing the fe-point correlation function through a determinant of size 
fc X A: of a single kernel. Our goal is the computation of this kernel. 
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3 Solution of the Half-Degenerate Case 

In this section we present the solution of the simpler ensemble (12.81) . i.e. in the half-degenerate case 
with Sn = a a^n- In Section [4] we will discuss the more general non-degenerate case. 

Let us first restate the starting point for our problem, the joint probability density (j2.8p : 

- , Xn) = CN^Ni^NB det[y5fc(Ai)|i<fc^K7v] • (3.1) 

Here we have simply used the second form of the Vandermonde determinant (12.9p and moved the 
Laguerre weight factors jjj gq^ (|2.8p into the rows of the second determinant comprising iFi, 

yielding the one-point weight (pj defined in eq. (I2.15p . 

The building blocks in eq. (|3.ip can be simplified further. In appendix where an alternative 
derivation for the joint density eq. (12.8p is presented, an identity for the confluent hypergeometric 
function iFi is derived. With the help of this identity the one-point weight pj from eq. (I2.15D can be 
expressed for these parameter values in terms of elementary functions: 


, (Na + Nb- N)\{Nb -N + k)'. 


XNa 


-1-k 


Na-1 

nW = E (-i)‘ j 

Nb-N 


+ exp[-cj5]A] (-1) 


-N-k 


k=0 


,^{NA + NB-N)\{NA-l + k)\ A^s 

kl(NB -N-k)\{NA-1)1 - as])^A+k ' 


(3.2) 


This relation also follows from the handbook on special functions |32j as explained in appendix O It 
has to be compared with the expression as a sum over incomplete Gamma-functions given in m, that 
simplifies to eq. (13.2p only for Nb = N. More importantly, in appendix O the normalisation constant 
is found in a closed form in comparison to the determinant of 2 F 1 in eq- (I2.10h : 


^N,Na,Nb 


a 


-NaN 


N 

n 

k=l 


a 


.N-Nb-1 

Bk 


-A^!AAr({(Tsj}) 


m-i 

n 

^ 1=0 


{Nb - N)\ 


{Nb-N+ 1)\{Na + Nb-N)\ 


(3.3) 


Because of the bi-orthogonal structure of eq. (13.ip it is well known, see e.g. [22], that the fe-point 
correlation functions defined in eq. (|2.13l) can be expressed in terms ol a, k x k determinant. 


Rk (Ai, ■ ■ ■, Afc) — det \KAi{Xii ^j)\i<i,j<k] • (3-4) 

For k = N we obtain an expression for the joint probability density which underlines the fact that 
this ensemble represents a determinantal point process. The kernel Kn{\,pl) in eq. (13.41) is in general 
given in terms of bi-orthogonal functions. 


N-l 

E 

j=0 


Kn{X,p) = '^'il^j{X)(t>j{k‘) , with 6jk= / dAV’j(A)^fc(A) . 


(3.5) 


Here 5jk is the Kronecker delta. The 4>j{X) are in the linear span of the set {pi{X),l = 0, ...,j} 
and the ipj{X) are polynomials of degree j and thus in the linear span of the monomials inside the 
Vandermonde determinant (|2.9p . 

In our case the set {^^/(A)} is special in the following sense. The functions defined in eq. (|2.15p are 
the same functions in A for all j and differ only through the argument {(J )^Thus in the following 
we construct polynomials orthogonal to the same function for different parameters {a^ — Ugj). This 
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will be done by expressing them through the expectation value of a single characteristic polynomial of 
H = AA^ + BB^ in subsection 13.21 with the help of the supersymmetry method |26H28] . The map to 
this expectation value is done by using the notion of the inverse Gram matrix in the next subsection l3.ll 
After computing also the expectation value of an inverse characteristic polynomial in subsection 13.31 
we conclude this section by discussing the result for the spectral density in subsection 13.41 There it is 
illustrated and compared with Monte Carlo simulations for two examples. 

3.1 The A:-point density correlation functions and their kernel 

In this subsection we use an alternative form of the kernel A'jv(A, /r) compared to eq. (13.51) . which follows 
from [25]. There the following form for the kernel was derived for general bi-orthogonal ensembles: 

N 

KN{x,y) = V’jiv) ■ (3-6) 

ti=i 


It contains the inverse of the Gram matrix defined as 


9ij 




dAA* ^y?j(A) = + i) 2 T 1 (m + 1 


NA,m + i;m + l; (cr^^ 



(3.7) 


where the second equality follows from an elementary integral. In [T3| the normalisation constant in the 
form (I2.10p was expressed through the determinant of this Gram matrix, = A^ldetp^], 

using the standard Andreief formula m, see eq. (|B.ip with k = I = 0. 

Our task here is to map the problem of inverting the Gram matrix, see eq. (13.61) . to the evaluation 
of the expectation value of a single characteristic polynomial, yielding orthogonal polynomials in the 
sense of eq. (13.51) . 

It is instructive to rederive the determinantal expression for the A:-point correlation function ()3.4[) . 
as it will help us to evaluate the kernel. Following the generalised Andreief formula, see |36j and 
our appendix [B] the fc-point density correlation functions can be written as an (N + k) x {N + k) 
determinant (|2.I4p . Leaning on this determinantal expression we use the following identity for block 
determinants. 


det 


o 

c 


d 

b 


det [a] det [6 — ca , 


(3.8) 


where o, b, c and d are matrices with a invertible. Identifying b = Okxk as the k x k matrix with zero 
in each matrix entry in the identity (13.8p . we find the form of the kernel (13.6p as a consequence. 
Moreover, using the same line of computation backwards we obtain for the kernel: 


KNix,y) = -C 


N,Na,Ni 


N\ det 


= NO 




N,Na,Nb 


N i-c 

ni 


i;rdAA-vaA)|-;;:^ 




0 


dXi det 


A^' 


l2<i<7V 


det 


v^jiy)\ 

(Ai) I 


^<3<N 

^<j<N 

2<i<N 


(3.9) 


Expanding the first line with the help of the identity (13.8p . it becomes immediate that the right 
hand side is indeed the kernel ()3.6I) . From the second line in eq. ()3.9p . the well-known identification 
Kn{x,x) = Ri{x) immediately follows. 

Starting from the second line of eq. dSSI), we can use a simple Laplace expansion of the second 
determinant with respect to the first row containing the ^Pj{y)- The x-dependence in the Vandermonde 
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determinant can be rewritten as a smaller Vandermonde determinant in A' = diag(A 2 , • • •, Atv), times 
a characteristic polynomial with A' as a matrix and x as its variable. Thus, we obtain 


KN{x,y) 


N 


N 


j\T 

'-^N,Na,Nb 


n/ Vj(y)det (^fc(Az)l 


2<1<N 


1=2 ■ 


J=1 

N 


X^N-l{{^ 2 , • • • , Atv}) ]^(Afc — x) 


k=2 


Gj( det[xl7v-i ^ n-i,Na,Nb-i 
i=i 


(3.10) 


where the functions (pj are still given by (I2.15P or alternatively (13.211 and the constants in the sum are 
defined by 


Gj = (-1)^+^V 


C 


SaVs 

N,Na,Ne 


iNB-N)\ 


y/ yf 

'■^N-1,Na,Nb-1 


{Nb - l)!(iVA + Nb- N)\^Na^Nb-n+ 1 lA. 


n 


(3.11) 


¥j 


In the second equality of eq. (j3.10p we have used that up to normalisation the integrals under the sum 
correspond to the expectation value of a single characteristic polynomial of an (V—1) x (N— 1) random 
matrix H' with correlation matrices = ga^n-i and - = diag((Tsi,..., crsj+i,... ^gbn) 

where aBj is omitted. Note that also Nb gets shifted to Nb — 1 to guarantee that the parameters in 
the one-point weights ipj remain the same. In particular we use the fact that N always appears in 
the combination Nb — N. The computation of the kernel is thus reduced to the computation of the 
expectation value of a characteristic polynomial, with the expectation value defined in eq. (|2.7I1 . 

The kernel (j3.inp together with eq. (j3.4ll is our first main result. In order to state the complete an¬ 
swer for the kernel we already give the result for the expectation value of the characteristic polynomial 
here: 


(det[xlAr_i H 


= NaKNb - 1)! / 

J-fi 



dz2 


{x- ZiaA- Z 2 (JBk) ■ 
(3.12) 


The derivation of this quantity will be done in the next subsection 13.21 using supersymmetric methods 

In the form (j3.10p the kernel is amenable to take the large-V limit, replacing the sum by an integral 
and the summands given by the constants (13.lip , the double contour integral (I3.12p and the functions 
(I2.15|) by their asymptotic values. The latter is known for general parameter values including large 
arguments and indices, see e.g. [32]. Hence the large-V limit of all correlation functions (13.41) can be 
derived in this way. 

Before coming to the expectation value of a single characteristic polynomial let us interpret the 
alternative representation of the kernel (I3.inh in comparison to eq. (13. 6 j) . Typically the latter is 
simplihed by making a change of basis of the following kind. Linear combinations within the linear 
span of the elements of the two determinants in eq. (13. ip are sought, such that in the new basis the 
Gram matrix ()3.7p becomes diagonal and thus easy to invert. This construction reduces the double 
sum (|3.6h to a single sum over functions that are bi-orthogonal. It is also known that one of the two 
sets of functions, the bi-orthogonal polynomials, are spanned by the monic powers A]~^ and are given 
by the expectation value of a single characteristic polynomial. The difference between the standard 
literature and our case is that here all polynomials (I3.12p are of the same order, namely of order 
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— 1, but nonetheless they build a basis. Usually one has for each order 0,1, 2,... one polynomial 
of this degree. The reason for this difference is due to the fact that all one-point weights ipj only 

differ by the parameter ^<7^^ — nothing else. Keeping this basis {‘fj}i<j<N in the second 

determinant of the joint probability density (13.ip . we have to get the same polynomial apart from an 
additional argument cbj, due to symmetry reasons. Hence, our result (j3.10l) is equivalent to build new 
polynomials of the same degree only from the linear span of the monomials while keeping the 

functions <Pj{y) untouched. Consequently the two sets of functions satisfy the following orthogonality 
relation, 

/*oo 

J dxP^l^{x)ipk{x) = G~^5jk . (3.13) 

Here we have introduced the polynomial of degree — 1 in monic normalisation that depends 

on all asi with i j, 


— ( det [xIat-i H ]')— G- ^{g 


yi/ y/ 


N 


-1 

i=l 




(3.14) 


This expression for the polynomials bi-orthogonal to the functions ipj resembles the Heine formula [M] 
for orthogonal polynomials. The normalisation constants Gj appear in this relation due to the diago- 
nalisation of the Gram matrix g. 

Equation (I3.13P can be easily cross checked by explicitly writing the expectation value in terms of 
the defining integral. The integration over all variables x, X 2 , ■ ■ ■ An times the Vandermonde deter¬ 
minant anti-symmetrises in all variables. Whenever j ^ k we will encounter the product (pk{\i)^k{x) 
in the Laplace expansion of the second determinant, which is symmetric in the two arguments A; and 
X and thus vanishes. 


3.2 Expectation value of a single characteristic polynomial 

We now compute the expectation value of the characteristic polynomial of degree V, 


Pn{x) = ( det 


x1n-WW^ 


\pA,'^B 
/N,NA,Nb • 


(3.15) 


From this expression the polynomial P^^_i (x) needed for the kernel (|3.10p simply follows by reducing 
N —>■ N — 1, Nb —)• Nb — 1 and omitting the eigenvalue aBj of We directly consider the 

non-degenerate case, T,a / aA^N-, as it is not more complicated than in the half-degenerated case. 
Moreover, the computation in this subsection sketches already the main ideas to be applied in sectionHJ 
where the generating function for the fc-point density correlation function is computed in the general 
case. This generating function is given by expectation values of k ratios of characteristic polynomials. 

Our derivation uses the duality between VUVUl and in the first step, then we express the 

determinant through Gaussian integrals of fermionic nature and perform the Gaussian average over 
W. In a final step we make use of the superbosonisation formula [381442| to reduce the remaining 
number of integrals to a minimum. 

The duality between WW^ and hU^VU implies that we can write eq. ()3.15p equivalently as 


Pn{x) = x^-^'^(det 




WW ) 


^A,^B 
N,Na,Nb ’ 


(3.16) 


as the matrix WW^ of dimension N and the matrix W^W of dimension Nw = Na + Nb have the 
same N non-vanishing eigenvalues. In the case that WW^ has full rank N then W^W has N\y — N 
zero eigenvalues. 
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It is well known that the determinant of a matrix can be expressed through an integral over 
Grassmann variables, also called Berezin integral. Let us introduce a set of Nw complex (anti¬ 
commuting) Grassmann variables Vi, with the following convention for complex conjugation, 


{vi,Vj} = 0, {vi,v*} = 0, {v*)* = -Vi, '^i,j = l,...,Nw 


with {,} representing the anti-commutator. 
We arrange the Vj in two sets of vectors. 


V = 


Va 

Vb 


( 


, ^^4 = 


Vl 




and Vb = 


^a^a+i 




(3.17) 


(3.18) 


with scalar product 


Nw 


(3.19) 


2 = 1 

The standard flat integration measure for Grassmann variables can be defined for one particular 
Grassmann variable Vi by only two properties, with an obvious extension to multiple integrations: 


/ 


dviVi = 1 , 


J dvi = 0 , V i = 1 ,..., Nw and 


Nw 

[dV] = dv*dvi. 

i=l 


(3.20) 


Here the last definition is required because of a possible minus sign. With these dehnitions the 
characteristic polynomial (I3.16P can be written as 


Pn{x) = 


- ^N-Nw 


[dV] ^exp -V\xtN^ - W^W)V ^ 

x^-^'^detpAj-^^detpBj-^s J g 


Sa,Sb 

N,Na,Nb 


-xVW 


X det 


^-1 


ItVa + Hat <8) VaPa 
liv O VgGj 


1n(^VaVI 


1 -1 


A''B 

® Inb + tN<s> v^vl 


(3.21) 


In the second equality we have performed the average in W according to the correlated Gaussian 
averages m- Note that the matrices VaV^, ^b^a matrices composed as a dyadic tensor by 

taking the dyadic product of the two vectors, which are bosonic objects. Hence its matrix entries are 
commuting. 

In order to simplify the determinant in the integrand of the block-matrix 


V = 


a d 
c b 


^ I ® Ina + Itv ® V^vl 


In^V^V^ 


(g) Inb + liv <8) VbH/ 


B 


we apply the identity detP = det[a] det[ 6 ] det[l — b ^ca ^d]. The first two determinants are 
det [a] = det ® Ina + ^Va^a = detps^]”^^ det iN + VjyjiJ^A , 

det[ 6 ] = det (g) Itvs + Itv G) Vglg = detps^]”^® det Iat-|-V^ l/gSe 


-1 


(3.22) 


(3.23) 


In the second step we pulled the matrices Ya and Yb out of the determinants and used a well-known 
duality for Grassmann variables det[lAr^ -|- 7 VaV)|] = (1 -|- 7 V^Va)“^ with 7 an arbitrary constant. 
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Note that vJVa and V^Vb are scalars as mentioned in eq. (13.191) . This fact allows to reduce the size 
of the determinants to N instead of NNa and NNb, respectively. Furthermore we may write in eq. 


d = Itv ® = (Ijv ® Va) (liv 0 Vl) , 

c = (g) VbVI = (Ijv ® Vb) (liv ® Vl) . 


(3.24) 


This becomes useful for the remaining determinant for which one can make the following expansion 

Tr(C)^' 


det[l — b ^ca ^d] = det[l — C]= exp 


-E 

L k=l 


k 


(3.25) 


The last factor in C is (Iat <8 ) V^) from the matrix d which will be moved to the first position of C, 
using the cyclicity of the trace at the expense of a minus sign. The minus sign is due to the fermionic 
(anti-commuting) nature of d and of b~^ca~^. This minus carries over to the inverse of the determinant 
det[l — C]~^, where C = db~^ca~^. We thus obtain 


det[l — b ^ca = det 


1 -1 


= det 


Iat — (lAf ® V^) ^Iat ® lAfs + Sb ® ® Vb) (3.26) 

X (Itv <8) Vl) (^Itv <8) Ina + ® VaVa) ^ (^a «) V^) 

tN - (liv + vIVb^bY' (1n + vIVa^aY^ V'iKiEA 


1 -1 


In the last step we have commuted the following two factors. Expanding the geometric series of the 
second factor into its von Neumann series we obtain 


(Itv <8) V^) [Itv ® ItVb + Eb ® V^V^ 


-1 


(liv®V^)J](-l)"(EB®VBV^)" 

oo 

J2i-vlvBJ:B)H^N^vl) 

k=0 

lAf + V^Vb'SbJ (lAf <8) Vb) ■, 


(3.27) 


and likewise for the two factors (lAf <8 * Va) and ^Iat <8) Ia^^ + Ea <8 * Vj^vJ^ as in eq. (j3.27p . We can 
now insert eqs. ()3.23l) and ()3.27l) into eq. (I3.2ip and obtain 


-1 


Pn{x) = detpA] ^'^detps] J i^V] e ^^^^det[a] ^det[6] ^det[l —6 ^co ^d] ^ 


= ^N-Nw 


[dV] e-^^^^det 


liv + V^>aEa + V^Vb^b 


(3.28) 


after cancelling all normalisation factors and combining the three determinants. This is the result for 
the expectation value of a characteristic polynomial valid for arbitrary covariance matrices Ea and 
Es. 

Moreover, in the case when these two matrices commute and thus can be simultaneously diag- 
onalised, which is in particular true for the half-degenerate case Ea = cta^n, the result (I3.28P can 
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simply be expressed in terms of the two sets of eigenvalues {cTAk} and {ask} of the two matrices T,a 
and T,b, 

r ^ 

PNix) = / [dV] n (^ + + vlVBCJBk) ■ (3.29) 

k=i ^ 

In this case the determinant reduces to a simple product. 

In the final step we apply the superbosonisation formula [381l42] . In the present case one can readily 
understand this formula via a Taylor expansion. For any entire function, f(z) = /k\, 

that only depends on the scalar product (j3.19p . / = f{V^V), its Grassmann integral can be efficiently 
evaluated. From eq. (13.201) it is clear that only the term containing all pairs v*Vi contributes. This 
term is contained only in the power , with multiplicity N\y\. Thus the Grassmann integral 

becomes in this case 

/ /• r r] 1 

d[V]f{VW) = j d[V]f^^^\0)l[vlv, = i-l)^'^Nwlf • (3-30) 

k=l ^ 


The integration contour 7 encloses the origin in positive direction. This expression is the superboson¬ 
isation formula in the one-dimensional fermionic case. The only complication from eq. (I3.29P is that 
the scalar products V^Vj^ and appear with different factors. Hence we apply eq. (|3.30p twice, 

for each of the two sets of variables Va and Vb- We arrive at 


Pn{x) = Na'.NbI <f 

J-fi 



dZ2 

2'Ki 


qZi+Z2 


Na+1 Nb+^ 

Zi Z 2 


det [xItv 


ziT^a - Z2P^b] ■ 


(3.31) 


Note that we rescaled the contour integrals by the variable —x to absorb the additional prefactor 
xN-Nw and the sign in the integrand. 

Equation (I3.3ip is the second main result of this section. The expectation value of a single char¬ 
acteristic polynomial reduces to this simple expression, valid for commuting correlation matrices Tia 
and 'Pb with = 0. From this the polynomial P^^_i{x) in eq. (j3.14p . trivially follows, by 

setting = ga^n-i j P'b — diag(cJBi,..., (Jbj-i, ctbj+i, ... -,(7bn) and by replacing N ^ N — 1 and 
Nb -a Nb — ^ in the average. Then this polynomial reads 


p^£,{x) = NA\{NB-iy.f^ ^ 


dZ2 


3 ^ 1+22 


2Tri 


{x - ZiaA- Z 2 (TBk) 


(3.32) 


and together with eq. (|3.1UI) determines the kernel and thus all fe-point correlation functions of this 
model. 

We conclude this subsection with a consistency check. As we have pointed out earlier in eq. ( 12 . 121 ) 
for equal correlation matrices our setting reduces to a single correlated Wishart-Laguerre ensemble 
of matrix dimensions N x N\y. In the completely degenerate case Pa = Pb = (7a^n it is well 
known that the expectation value of a characteristic polynomial is simply given by the generalised 
Laguerre polynomial orthogonal with respect to the weight function w{x) = x^^~^ ex.p[—x/a] in 
monic normalisation, cf. eq. (12.lip for the joint density of all eigenvalues. Our result agrees with this 
limiting case as follows. 

In eq. (|3.29p a single application of the superbosonisation formula (j3.30p suffices when setting all 
eigenvalues of the correlation matrices Pa and Sg to be equal, aAk = (XBk = < 7 . From there we obtain 
for the fully degenerate case ( 12 . 111 ) 

pfy9^x) = Nwlf^^^jj^{x-za)^. (3.33) 
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This has to be compared to the standard complex contour integral representation of the generalised 
Laguerre polynomial, see e.g. [32] . 


L“(x) 


I dz 

Jc 27ri (z — x" 


(-ir 


dv {xa — ti(T)"+“ 

27ri yri+l 


(3.34) 


In the second equation we have simply shifted and rescaled the contour C, formerly enclosing the point 
X and not the origin, to a contour 7 around the origin (and not including v = x). After an appropriate 
rescaling we thus have 





Nw 

Nw 


Nl 


^ ^ ^^{Nw-k)\T{k + N-Nw + l)k\ 


X\k 

a/ 


N 


i-afmY, 


Nwl 


1=0 


{N - l)W.{Nw - N + ly. \ a 


i-afNlL^'^-^ix/a). 


(3.35) 


In the hrst step we have used the explicit representation of the generalised Laguerre polynomial. Due 
to the fact that N — Nw = N — Na — Nb < 0, the Gamma-function in the denominator truncates the 
sum from below up to A: = Nw — N. A shift I = k — Nw + N leads to the desired form, the generalised 
Laguerre polynomial of degree N in monic normalisation. It is orthogonal with respect to the weight 
function as desired for the limiting ensemble (| 2 .Iip with a = aA- 


3.3 Expectation value of an inverse characteristic polynomial 

Let us compute the expectation value of an inverse characteristic polynomial 

QN{y) = (det[yl,v - 


N,Na,Nb 


= y-^ (det[l^^ - 


N,Na,Ne 


(3.36) 


Here we have to choose Im(y) 7 ^ 0 in order to regularise the expression. This quantity is related to 
the Cauchy-transform of the one-point weights ipj. This can be seen by combining the Vandermonde 
determinant A 7 v({Aj}) of the joint probability density (13.11) with the inverse determinant in eq. ()3.36p . 
i.e. 


^Af({Aj}) 

nf=i( 2 /-A, 


= det 


5 . A;-I|l<fc<v -1 
\i<j<N 


1 


y-^j 


^<j<N 


(3.37) 


see |3B]. This determinant can be expanded in the last row yielding N terms in eq. (13.11) . All N terms 
give the same integral such that 


N 


QAy )=n / 

j=i Jo 


dX 


2Xn-i{Xi, ..., Ajv_i) 

y - Xn 


det 


V’l[Xk)\i<i<j\f 

^li^N)\l<l<N 


(3.38) 


Likewise we expand the second determinant in ipi {Xn)- The remaining integrals yield the constants 
1/C^^^^n^^^Nb-i- Then we have 


./V poo J /* 

QAy) = = I 

• 1 ^ J—QQ y j_ 

'-'N-1,Na,Nb-1 


dx 


N 


-ooy-x 




(3.39) 
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Here, we used the definition (j.S.llj) for the constants Gj. Hence QN{y) is the Cauchy transform of a 
weighted average of the functions (pj (x). Note that in the standard setting of bi-orthogonal polynomials 
QN{y) would be given instead by the Cauchy transform of a single orthogonal polynomial. 

Let us come to the calculation of the average (I3.36p . In this case the inverse determinant can be 
expressed as a Gaussian integral over two complex vectors Ua and Ub of ordinary bosonic variables, 
in analogy to eq. (13.181) . Apart from signs due to the fermionic nature of V and the bosonic one of 
U the computation is very similar to the one presented in the previous subsection. In contrast to the 
fermionic case we will encounter poles, which is why we have to keep track of the regulating imaginary 
part Im(y) / 0 . 

After integrating over the Gaussian matrices A and B, using the duality (j3.36p and identities for 
determinants we arrive at 


QN{y) 




I [dU] 


e-^'^det 


y^N 


u\Ua^a - U^Ub^b 


(3.40) 


with [dU] being the flat measure on . This is valid for arbitrary covariance matrices T,a and In 
the case of commuting covariance matrices, = 0 , and in particular in the half-degenerate case 

the determinant in eq. (j3.40l) can be diagonalised, and the integral simplifies further. The bosonisation 
formula for the present case reads 


/d[t/l/(t/tt/) . (3.41) 

0 

This follows from going over to polar coordinates in 2Nw real dimensions for U, with the scalar 
product U^U = s being the squared norm of the complex vector U. Applied to both Ua and to Ub 
we find the hnal expression valid for commuting covariance matrices, compared to eq. (|3.31l) : 

OO OO 

qr^w r r 

QN{y) = _ i)!(Arg - 1)1 J J _ 52 ^ 5 ]"^. (3.42) 

0 0 

The correct monic normalisation can be easily checked, by observing that lim|y|^oo y'^Qjv(y) = Ij as 
it is required from the definition (13.361) . 


3.4 The spectral density 

In this subsection we will apply the previous results to the spectral density and give an alternative 
integral representation for the density and kernel. They are illustrated by numerical simulations below. 

For this purpose we combine the result for the orthogonal polynomial (13.321) with the expression 
for the kernel (j3.inj) including (pj, see eq. (j2.15p . In particular we exchange the finite sum with the 
integral and have 


KN{x,y) 


^yNA+NB-N^-y/cTA 



dZ2 e^l+^2 


N 

iFi (m + 1 
i=i 


NA-,m + 1 ; (o-^^ 



(3.43) 


I 


Nb-N+I 

^Bj 


n 


{x - ZiaA- Z 2 CFBl) 
{aBj - ctbi) 


with the constant 

NaI{Nb-N)\ 

"" {Na + Nb-N)\' 


(3.44) 
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Figure 1: Comparison of the analytical result (|3.46p for the spectral density (red curves) with 
Monte Carlo simulations (histograms=shaded area, 10® matrices drawn from the ensemble (| 2 . 1 l) i. 
We employed the parameters N = 9 with the fixed covariance matrices = Ig and T,b = 
diag(0.02,0.20,0.30,1.50,2.01,2.25,2.27,4.05,4.13) and the time length of the epochs {Na,Nb) = 
(35,40) (left plot) and {Na,Nb) = (350,400) (right plot). For longer epochs Na and Nb the peaks 
overlap less compared to the shorter epochs. Hence they shift more and more towards the determinis¬ 
tic positions given by equation (j3.47l) . The remaining deviations from those positions result from the 
level repulsion amongst the individual eigenvalue distributions originating from the joint density (j 2 . 8 p 
which are still visible. 


From eq. (|3.43p it is obvious that the kernel is expressible in terms of matrix invariants of T,b- Therefore 
we have to extend the product in eq. (j3.43p by the missing terms in cr^j. This can be achieved by 
introducing a contour integral in an auxiliary variable as 

K^(X V) =cv^A+N,-N^-y/.^ I _1_ 

KN{x,y) cy f ^ ^ 27 rz det[z3lw - 


71 


72 


(3.45) 


73 


det[{x - ziaA)tN - Z2 ^b] ^ , , , ,,, , ^ \ 

--j-iFi (m + 1 - Na \m + 1; [a^ - z^jy) 

X - ZiaA - 22^3 


which replaces the sum. Here the contour integral over 73 has to be specified in the following way. 
The contour 73 only encircles the positive eigenvalues of but excludes the pole Z 3 = Z 2 /{x — z\aj^. 
Since z\ and Z 2 lie on the contours that encircle the origin and no other pole, we can choose the radii 
of these contours equal to [x — e)l(yA and eminj=i^,.._w 


a 


Bj 


/2 


respectively, with x > e > 0. Then 


a 


-1 


Bj 


In this way we have 


the contour 73 encircles the positive real axis starting from minj=i^...^jv 
the desired realisation of the contours given in the representation (j3.45p . Thereby the poles from the 
determinant det[^; 3 ljv — yield the only contributions to the integral over 73 . 

The expression (I3.45|) is suitable for an asymptotic analysis at large matrix dimensions, as an 
alternative to the discussion of eq. (|3.10p at the end of section 13.11 Here one has to deform the 
contours suitably such that the saddle-point analysis can be performed. We will not pursue this 
further as the aim of the present work is the derivation of exact results at finite matrix size. 

Following eq. (j3.45p the spectral density is given by 


Ri{x) = Kj\f{x,x) 


^^Na+Nb-N^-x/cta 


71 


dzi 

2 Txi 


72 


dZ 2 

2 Txi 


dzs 


^Zl+Z 2 


73 


2vrz det[z3lw - ^ b ^' 

det[(x - 2 iCrA)lAr - 22 Sb] „ / , ^ , 1 / -1 ^ \ 

:-^1 iFi (m + 1- NA;m + l; (a/ - Z 3 )x) 

X - ZlUA - Z2Z^ 


(3.46) 
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Since all integrals can be evaluated by the residue theorem at the specific poles, the simplest way to 
compute this three-fold integral numerically is via its series expansion. In Fig. [T]we show the spectral 
density for hnite N for two generically chosen sets of eigenvalues of Here we set cta = ^, since a a 
only rescales the spectrum. The density is compared to numerical simulations, with details given in 
Fig. [H and we find an excellent agreement. 

Let us try to understand some qualitative features of the density shown in Fig. [TJ The position 
of the peaks of the spectrum can be easily estimated in the regime Na,Nb S> N. Indeed when 
scaling Na = nN'^, Nb = and x = nx' with N,N'j^,N'q,x fixed, in the limit n —)• oo we 

can perform a saddle-point approximation of the matrix representation of the ensemble in eq. (ET}, 
yielding H = uN'^^a + nN'^T^B- Then the limiting spectral density simplifies to 

lim nRiinx') = Tf <5 (x'liv - N'^Y^a - N'bYb) ■ (3.47) 

This limit also holds in the general case where both matrices Ya and Yb are non-degenerate. At finite 
but large Na^ Nb ^ N the peaks are broadened. When the distributions of the individual eigenvalues 
of H overlap the eigenvalues repel each other and shift away from the deterministic positions (I3.47jl . 
This can be nicely seen in Fig. [H where we considered two examples with (Na,Nb) = (35,40) and 
with (Na,Nb) = (350,400), keeping N = 9 fixed. The sharpening of the peaks for larger values 
of Na and Nb is obvious. Nonetheless the distributions of the individual eigenvalues still overlap. 
In a Gaussian approximation the spacing between neighbouring peaks scales as Na + Nb while the 
width of the peaks scales as \/Na -|- Nb- Thus a factor of ten in Na and Nb yields a factor of about 
l/VTO ~ 1/3, which explains the shape modification form the left to the right in Fig. [TJ 


4 Solution of the Non-Degenerate Case 


In this section we will compute the generating function for the fc-point density correlation function in 
the general case of non-degenerate Ya Yb- In view of the distribution of matrix elements of H being 
given by a hypergeometric function of matrix argument, see eq. ra. we have to choose a different 
strategy, due to the absence of any bi-orthogonal structure. The techniques we will use instead are 
supersymmetry and superbosonisation, generalising the computations from subsections 13.31 and 13.21 
Let us hrst redehne the fc-point density correlation functions, 

Ik \ 

i?fc(Ai,...,Afc) =/fjTr 5(A,-1^-77) \ . (4.1) 

' 7=1 / N,Na,Nb 

The matrix delta-functions in this expectation value can be generated by an appropriate differentiation 
of the following generating function for the fc-point density correlation function. 


Z,\p{X) 


/ n^=idet[x,lAr-H] \^^’^" 


(4.2) 


where we denote X = diag(yi,...,xi,...,Xp), with lm{yi) / 0 for all I = 1,... ,q. The generating 
function Zq^p{X) defined in eq. ()4.2p is the central object of this section. 

As an important example let us consider the spectral density i?i(Ai) with k = 1. In order to derive 
it let us define the averaged Green’s function or resolvent as 


VFi(y) = 



1 

ylN - H 


N,Na,Nb 


(4.3) 
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It is well known that the spectral density can be obtained from the resolvent by taking its imaginary 
part: 

Ri{y) = - lim lm{Wi{y)). (4.4) 

Vr Im(j/)—>-0+ 

Note that the density i?i(Ai) is normalised to N. From eq. (14.21) at p = q= l = k the resolvent Wi{y) 
results from the generating function by simple differentiation: 

d,Z,^,{X)l^^ = W^iy). (4.5) 

In the same way the product of all k Dirac delta-functions leading to the A:-point density correlation 
functions Rk{Xi, ■ ■ ■, Xk) in eq. (14.11) can be generated from Z^k{X) by successive differentiation leading 
to the A:-point resolvent, and by taking imaginary parts to zero from above as in eq. (|4.4p . 

It is well known [T] that for k > 1 the /c-point density correlation functions defined in eq. M 
differ from those defined in eq. ( 12 . 1 ,ip by so-called contact terms, e.g. for k = 2 

-^ 2 (^ 1 , A 2 ) = i? 2 (Ai, A 2 ) -k 5(Ai — A 2 )i?i(Ai) , (4.6) 

where the last term originates from coinciding arguments ot k = 2 Dirac delta-functions in eq. (HU). 

We would like to add that for <7 = 0 in eq. (|4.2I) the quantity Z^^p can be used in combination 
with the replica method to compute the generating function of the ergodic capacity, the average of 
ln[det[l -k H]]. As mentioned in the introduction this quantity is of central interest in applications to 
telecommunications in the setup of MIMO. 

The strategy to compute all ratios of characteristic polynomials in the generalised partition function 
Zq\p{X) will be similar to the computation in subsection [3]2l We rewrite the determinants as Gaussian 
integrals, now over commuting and anti-commuting variables, then integrate over W, apply the duality 
between ordinary and supermatrix spaces as in eq. (j.‘l.,S 6 p and in the end employ the corresponding 
superbosonisation formula |381t42j . The details of this computation are carried out in appendix [Cj 
The final result for q < N is 

Zg\p{X) = Cn^Cns J dp{UA) I (Ua) Sdet^^ (Ub) (4.7) 

xSdet“^(lAr ® X — T.a®Ua — Z.b® Ub) , 
where the constants are given by 


q-l 


c„=n 


p-1 

n 


{n- q + ly. 


-I--I- vrTn — o -k /)! 

1=0 ^ ^ ' 1=1 


7r‘ 


(4.8) 


The definition of the supermatrices Ua, Ub S Herm_|_(g|p), the supertrace Str(...) and the superdeter¬ 
minant Sdet(...) are recalled in appendix 1C. 11 The Haar measure is given in eq. (IC.lSp . 

As a check for Tia = = S, eq. (|4.7I) simplifies to the well-known supersymmetric result |12lll9ll43j 

of a single correlated Wishart-Laguerre ensemble, 

^q|p(^)|E^,EB=E = Cn^+Ns I dfi{U)e-^^^^Sdet^^+^^{U)Sdet-\lN^X (4.9) 

This result can be readily deduced from eq. (14.7p via the two substitutions Ub = U\^'^U'U\^^ and 
Ua = (Itv + f7')“^/^f7(lAr -k U')~^/'^. Here we have used the group invariance of the Haar measure. 

In the case of the spectral density we need the case {q\p) = (1|1)- Then the Haar measure is simply 
dp{U) = {27ri)~^[dU], with [dU] the flat measure. After a brute force expansion in the Grassmann 
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variables and the integration over the two angles the generating function reduces to the one shown 
in eq. (ID.281) . For the calculation we refer to appendix [Dj The derivative with respect lo x aX x = y 
yields the spectral density, 


My) = 


t:{Na-1)\{Nb-1)\ 

d ( d 

d 


(4.10) 


X lim Im 1 1 + 

Im 1 /—>0+ 


sasb 


+ — 1 + 


SA 


X det 


1 + 


d 


duA 


1 + 


d 

dfiB 


dyBj 
1 


92 


+ 


+ 

SB 


1 + 


poo 

dsA / det[Fs] 

Jo 


d 


dy-Bj dzAdz\ 


dfiAJ SzbBz*^ dzAdz*j.dzBdz 


B 


ylN - yA^A - ys'^B + [ ^A^A, Zb^B ] Fg 


Tr ( yl/v - yA^A - yB^B + [ zaT^a, zbT.b ] Fg 


■ 1 

i 


1 — 

1 — 

w 

0 

i 

1 


z*b'Zib 

■ 1 

i 


1- 

w 

1_ 

i 

1 


Z^ZjB 


-1 




with 


Fs = {y^N - sa'F‘A - sbF>b) ^ 


(4.11) 


The derivatives in yA/B za/b encode the integrals we carried out and can be easily numerically 
evaluated via a series expansion. The two-fold integral over and s_b is the numerically non-trivial 
part, since it does not factorise. The limit of vanishing imaginary part, Im(y) —)• O"*", reduces one 
of these remaining integrals to a Dirac delta-function. Thus we effectively end up with a sum of 
one-dimensional integrals, which can be evaluated numerically. A further and a more efficient way to 
evaluate the latter expression numerically is to adjust the imaginary increment. It can be chosen small 
enough for a result independent of the imaginary increment, while it is big enough for a numerically 
stable evaluation. 

The expression for the spectral density simplifies a lot in the case when Na, Nb and N become 
large. We will not perform a rigorous asymptotic analysis here but give a sketch of the expected 
answer for the limiting global or macroscopic density. Our approach complements the ideas for taking 
the large-limit sketched in section [3l For this purpose we follow the ideas of refs. |44p45j and start 
with the resolvent in the supersymmetric formulation 


VFi(y) = 


1 


Zi\i{yti\i 
xStr 


dyiUA) / dy{UB) 


e-Str^A-Strt/Bgjg^jVA Sdet^^ {Ub) 
Sdet(?/lAr (g) li|i - Tia®Ua-Fib® Ub) 


(4.12) 


{y\N — A® uA —F.b®Ub) ^ f q 


The normalisation can be indeed chosen as 1/Zyi[y^yi) which is certainly equal to 1, with the 
advantage that we see which terms in the saddle-point approximation cancel. The saddle-point ap¬ 
proximation starts with the Lagrangian 


C{Ua, Ub) = Str {Ua + Ub) — Aj^Str IiiI/a - IVeStr InC/e -I- Strln(?/lAr (g) li|i -'Za®Ua-'Ub®Ub), 

(4.13) 

where all superdeterminants have been raised to the exponent. From the first derivatives of the 
Lagrangian the saddle-point conditions, [7®) = 0 and dugF{U^\UQ^) = 0, yield 




- Tr 


N 


(yljv ® li|i - ® U^a'^ - Sb (g) U''s’)-^^A 




= 0 , 


(4.14) 
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and the second equation results from interchanging A and B. The notation TV n indicates the partial 
trace in the ordinary N x N matrix space, such that the result of this trace operation is an (1|1) x 
(1|1) supermatrix. The remaining trace over this supermatrix is denoted by Str(i|i). Assuming the 

uniqueness of the contributing saddle point {U^\ we have for the spectral density 


Ri{y) = 


—Im Str 

TT 

—Im Str, 

TT 


{ylN ® li|i - Sa ® - Sb ® UP)-^ 


In 

0 


(i|i) 


Tr 


N 


{ylN ® li|i - Sa ® - Sb ® UP)-^ 


1 0 
0 0 


(4.15) 


Note that we have already taken the limit Im(?/) —)• 0^ here which is possible due to the regularity 
of the saddle-point solution. To see that no Grassmann variables are involved in this expression 
(otherwise it would be inconsistent) we multiply eq. (14.141) with and add it to the corresponding 
equation for This leads to 


Ua + up - {Na + Nb- iV)li|i - yTr JV [(yl^v ® Iqi - Sa ® - Sb ® uP)-^] = 0. (4.16) 


We can insert this equation into eq. (I4.15P and arrive at the simpler expression 


Riiy) 


—Im Str 
vr 


(i|i) 


l(A°’ + A°’) 



(4.17) 


Hence only the diagonal elements of the saddle point {uP, uP) are involved. 

Analysing the saddle-point equation (j4.14p in more detail requires some effort. However, as we 
have already seen in the exact expression (j4.10p for the spectral density, the Grassmann variables 
only contribute to the integrand as a polynomial of a very small order. Therefore they do not enter 
the saddle-point analysis. Thanks to the projection matrix in eq. (I4.17p only the upper boson-boson 
entries of the matrices uP and uP enter, denoted by qP and qP, respectively. They are determined 
by the two coupled saddle-point equations for these entries, eq. (j4.14p and its counterpart: 


1 - — - Tr 
1-^-Tr 

yP 


{ylN - qP^A 
{ylN - qP^A 


qPj:B)-^^A 

qPj:B)-^^B 


0 , 

0 . 


(4.18) 


The two unknown complex functions qP {y) and qP {y) depend on the variable y. The asymptotic of 
the spectral density thus simplifies to the final result 


Ri{y) = Q(gA^( 2 /) + qPiy))^ ■ (4.19) 

This expression is reminiscent to the one presented in ref. |44] for the spectral density of a single 
correlated Wishart matrix. However, the evaluation of the corresponding saddle-point equation in |44j 
compared to those in eq. (I4.18P was much simpler because of two points. 

First, one is looking for only one function q{y) in the case of a single correlated Wishart matrix, 
satisfying a single and not two coupled equations. In our case there are {N -|- 1)^ solutions including 
multiplicity due to Bezout’s theorem, underlining the complexity of the problem. These solutions 
have to come either in complex conjugate pairs or are real because we can complex conjugate both 
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equations. Due to the imaginary increment of y only one half of the solutions can be reached with the 
integration of sa and sb, see the exact result ()4.10p . 

Second, only one empirical correlation matrix S is involved for a single Wishart matrix [44]. Thus 
we can diagonalise this matrix, which is not possible in the most general case for the present random 
matrix model with two empirical correlation matrices, ^ '^b satisfying [S^, 7^ 0. Nonetheless, 

the range of applicability of the strikingly simple result (I4.19P should be as good as for the simpler 
model studied in ref. jUj. We close the discussion at this point because the main goal of the present 
work is the exact finite-solution and not the detailed asymptotic analysis. 

Let us make a few hnal remarks on generalisations of the result (021). One can easily extend 
this to real (/3 = 1) and quaternion (/3 = 4) matrices H. The symmetries of the supermatrices 
change accordingly, where the positive definite Hermitian boson-boson blocks of Ua and Ub become 
either symmetric or self-dual and the unitary fermion-fermion blocks become self-dual or symmetric, 
respectively, see |41H43l[46] . Moreover, these supermatrices become twice as large for both (3 = 1,4. 
Another point, which only plays a role in the real case, is the exponent of the superdeterminants (14.7|) . 
For (5 = 1 one has to take the square roots of all three superdeterminants, while for /3 = 4 the 
exponents remain the same. 

Moreover, with the focus on a generalisation of the generating function (j4.7|) one can replace the 
Gaussian weights for A and B in eq. (12.ip by other weights from invariant ensembles like the Jacobi 
or the Cauchy-Lorentz ensemble cf. [471449| for recent related works. A simple way to calculate the 
according expression is shown in the works |43ll46j and is called projection formula. It is a shortcut 
of the map to superspace and does not only show that the weights in superspace (here the terms 
Qa{Ua)Qb{Ub) oc exp[—Strf/yi — Strf/fi]) have to be replaced, but also what these weights Qa and 
Qb are as functionals of the probability distributions Va and Vb in ordinary space. 

For another generalisation, more than two epochs can be introduced, namely the matrix H = 
AiA^^ + ... + with correlation matrices YIai , • • •, ^At and weights Va\ j • • • j Bat for arbitrary T 

can be considered. Then the result (14.71) can be extended to the product of T integrals, 

T 

= [t\CN^ MU A ,) Qa, {UA,)SdetM {Ua, ) (4.20) 

i=i 

X Sdef^ ( Iat (g) A - ^ Saj, G) Ua^ 

\ k=l 

As mentioned above, one could also choose weights Va, different from eq. m, e.g. by choosing 
Gaussian ensembles of Hermitian matrices or ensembles of Jacobian or Gauchy-Lorentzian type. Then 
the first line of ()4.2np changes accordingly in the supersymmetric weights Qa, , depending on the initial 
weights Va, ■ In the second line of eq. (I4.20p the substructure of the generalisation to a matrix H as 
the sum of T correlation matrices is clear. It is astonishing that our result for the generating function 
for H as a sum of T = 2 Wishart correlation matrices (j4.7|) enjoys a generalisation to much more 
complicated ensembles. 

5 Discussion and Outlook 

In the present work we determined the spectral statistics of all /c-point eigenvalue density correlation 
functions for the sum H of two correlated Wishart matrices. In previous work by Kumar it was 
shown that the joint density of the eigenvalues of H is given by a hyper geometric function of matrix 
argument. It was also shown in his work that in the half-degenerate case, when one of the correlation 
matrices is proportional to the identity, the joint density becomes a bi-orthogonal ensemble. It is given 
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by the Laguerre weight times a Vandermonde determinant and a determinant of ordinary confluent 
hypergeometric functions of Kummer type iFi. We first solved this half-degenerate case by determin¬ 
ing its kernel of bi-orthogonal functions explicitly, exploiting the bi-orthogonal structure. Remarkably, 
the kernel could be constructed by computing the expectation value of a single characteristic poly¬ 
nomial, providing the polynomials orthogonal to iFi. This expectation value was computed using 
supersymmetry and bosonisation in the general non-degenerate case. As a byproduct we computed 
the expectation value of the inverse of a single characteristic polynomial in this general case too, as 
well as a more compact form of the normalisation constant of the joint density in the half-degenerate 
case. The solution for the A:-point correlation function is then given in terms of a determinant of size 
k of this kernel. In contrast the previous representation derived by Kumar was given in terms of a 
determinant of size N + k, containing the moments of iFi. Our result thus offers the possibility to take 
the large-V limit and to study questions of universality by making an asymptotic expansion of the 
kernel. In addition to the sum over bi-orthogonal functions we gave an alternative three-fold complex 
contour integral representation of our kernel. Both representations are amenable to take the large-V 
limit. We illustrated our flnite-V result for the spectral density with Monte-Carlo simulations. 

In the more difficult case where both correlation matrices are non-degenerate no bi-orthogonal 
structure is available. Therefore we applied the standard supersymmetric method and computed the 
generating functions for the fc-point resolvents. They are given by the expectation values of k ratios of 
characteristic polynomials that we computed using the superbosonisation formula. They are given by 
two integrals over symmetric supermanifolds of fc-dependent dimension. As an example we spelled out 
the generating function for the spectral density by choosing an explicit parametrisation for k = 1. It 
can be written as an integral over two real positive variables and two angles, containing determinants 
and traces of rational functions of the correlation matrices and the four integration variables. The 
angular integrals could be solved, leading to derivatives of the order of the matrix dimensions. Based 
on these results we derived two coupled saddle-point equations in two variables in the large-A^ limit 
and expressed the limiting spectral density in terms of their solution. 

Several open questions are left for future work. This includes in particular the asymptotic analysis 
of the limiting kernel from the first part of our investigations, but also a detailed analysis of the 
saddle-point equations for the density in the more general non-degenerate case from the second part 
of our work. In this second part we also gave indications how our results could be generalised in 
several directions. This includes different symmetry classes with Wishart matrices built from real or 
quaternionic matrix elements. We also wrote down the general structure of the generating function 
when the matrix H consists of the sum of T > 2 correlated Wishart matrices in the non-degenerate 
case. This structure persists even when the Gaussian weights of the T >2 correlated Wishart matrices 
are replaced by more general weights of ensembles of Jacobi or Cauchy-Lorentzian type. Of course in 
this case the corresponding supersymmetric weights inside the integrals over the supermanifolds have 
to be replaced by their Jacobi or Cauchy-Lorentzian counterparts. 
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A Simplification of the Joint Probability Density 


In this appendix we present an aiternative derivation of the joint probabiiity density Pn, see eq. (12.8p . 
The benefit of this computation is twofoid. First, we obtain a more expiicit form of the normaiisation 
constant in eq. ()3.3I) given by a product of Gamma functions and the eigenvaiues of the non-triviai 
covariance matrix. This has to be compared to the determinant of a Gauss type hypergeometric 
function in eq. (|2.1Up . Second, we show that in the speciai case of parameters at hand the confluent 
hypergeometric functions inside the determinant of the joint probability density (12.8p can be expressed 
in terms of more elementary functions. 

To calculate the joint probability density of the matrix H = AA^ we introduce a test 

function f{H), which is a Schwartz function on the space Herm(A^) of Hermitian N x N matrices and 
satisfies the relation 

f{H) = for all H G Herm(A^) and U G U(A^). (A.l) 


Then for any function / with these properties we have that 

{f{AA^ + = I[dH] VH{H)f{H) = J [dA] iP^(Ai,..., A^)/(A). (A.2) 


In the sense of weak topology this is the definition of the joint probability density "Pat of eigenvalues 
A = diag(Ai,..., Xjsj) of the combined random matrix H = AA^ + BBK We make use of exactly this 
definition to derive an alternative expression for the joint probability density Pn in the half-degenerate 
case which is more explicit than the one derived in Ca¬ 
in a first step in deriving the joint probability density we introduce the Dirac delta-function 

N 

5{H - AA^ - BB^) = 5{Hjj - (AA^ + BB^)jj) (A.3) 

i=i 

X H 5iRe[Hij-{AA^+ BB%])5{lm[Hij-{AA^+ BB^)^j]). 


Those Dirac delta-functions can be written as a double Fourier transform yielding 

(/(AAt + Ppt))WB^^ (A.4) 

= S^^At] J[dB] exp[-Tr 

X [[dH]f{H) [[dK] exp[-tTV A:^] exp[iTr A:(P - AA^ - BB^)] . 


We also introduced a regularising Gaussian factor exp[—tTr with auxiliary parameter t, which we 
send to zero in the end. In this way all four sets of integrals over A, B, H and K are absolutely 
integrable and, thus, can be interchanged. The integral over H is absolutely integrable because of the 
test function /, which is also the reason for employing this kind of function. 

The first fraction in the constant in front of the integral (IA.4P is the one of the original probability 
weights ()2.1j) . The second fraction is the normalisation of the double Fourier transform such that the 
Dirac delta-functions are normalised to unity. 

In a second step we integrate over the matrices A and B. Hence eq. ()A.4p becomes 


(/(AAt + Hpt)) 


N,Na,Nb 


I [dH]f{H) I [dK]exp [-tTr K‘^ + iTr K H] 
xdef^^ + iK] def^s + iK] . (A.5) 
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To proceed further we have to simplify the ensemble since we want to diagonalise the matrices H = 
VAV^ and K = UQU^ and we need to integrate over the cosets V,U & U(A^)/[U^(1) x §(A^)]. The 
set §{N) is the permutation group of N elements and lifts the ordering of the eigenvalues Xj and Uj 
of H and K, respectively. 

At this point of the calculation we assume Y^a = ctaItv- Then we can also assume that T,b = 
diag(crBi,..., ctbtv) because the whole system is invariant under adjoint transformations —)• 

for all unitary matrices O £ \J{N). Diagonalising H and K and using the invariance of the normalised 
Haar measure of the coset U(A^)/[U'^(1) x S(iV)] under V —)• UV we have 


(f(AA< + 

det-^^SAdet-^^Y^B 


'•OO 

d^iU) exp[iTr C/nC/^A] 


J MA]A^({Aj})/(A) J 

d/i(l/)det-^s + iVnV^ 


(A.6) 


The additional constant is the squared volume of the coset U(A)/[U'^(1) x §(iV)]. The square is 
needed since we have one coset for U and one for V and the Haar measures are normalised. 

The two coset integrals are well-known. The first integral is the Harish-Chandra-Itzykson-Zuber 
integral [501IM] 


J dfi(U) exp[iTr C/H[/^A] 


‘i®t[exp[ia;aA6]|i<a,fe<Ar] 

A^({A,})A^({a;,}) 


where the normalisation is fixed to unity at H = 0. The second integral is also known |52j 


(A.7) 


/N-l 


dfi{V)det-^^ + iVnV^ = 11 " 


1\{Nb-N)1 


det 


{o'b\ ^\l<a,b<N 


\N-N„-li 


lA ii(^NB-N + l)\^ 




(A.8) 


Again the normalisation can be fixed by taking H = 0 which yields det ^^YIb- We insert these two 
integrals into eq. (jA.6p and find 


(/(AAt + 


(A.9) 


(27r)^[A!]2Ajv({aB,}) 


n 


{Nb-N)\ 


U ^Nb-N + 1 )\ 


J [dA]AN{{Xj})f{A) 


: J [dH]t 


-tTTn^..-NA r^-i 


det -|- iH] det 


iuJa I 


det 




N-Nn-1\ 


\ l<a,6<A^ 


To evaluate the integral over Q we apply Andreief’s integration formula [37], see eq. (jB.lj) for 
k = I = 0. Thereby we notice that the integral over H is absolutely integrable even at t = 0 because 
of Na, Nb > N^ such that we can set t = 0 from now on. Thus we end up with 


(/(AAt + 

det-^^SAdet^-^s-^Ss {Nb - N)\ 

NlANiWBj}) {Nb-N + 1)\{Na + Nb-N)\ 


(A.IO) 

I [dA]A^({A,})/(A) 


X det 


{Na + Nb-N)1 


dn exp[mAa] 

27r -|- iK)^B-N+i 


l<a,b<N 
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Now we are ready to identify the explicit form of the normalisation constant 


'-'N,Na,Nb 


det-^^SAdet^-^^-^SB 


/7V-1 

n 

V 1=0 


(Nb-N)] 


{NB-N + iy.iNA + NB-N)] 


(A.ll) 


of the joint probability density ■ ■ ■ ,^n) and the functions inside the determinant in eq. (lA.lOp 


^,-(A) = {Na + Nb 



dn exp[i«;A] 

27r [a^ + iK)^^{cry^] + iiPj^s-N+i ' 


(A.12) 


These functions are normalised such that lim^^-o = 1- They can be calculated via 

the residue theorem. Thereby we only get a contribution when A is positive. Then we can close the 
contour around the two poles and This yields two contributions, and we have explicitly 




( \ ^A — 1 1 

d \ exp[-c7j^U] 

d<yA) K'- 

_ d_\ exp[-ag^.A] 

9 ^ 3 ]} 


iVA-l 

exp[-a2^X] {-N 


k=0 


{Na + Nb- N)1{Nb -N + k)\ 
k\{NA - 1 - ky.{NB - N)] 





Nb-N 

+ exp[-a;j^]A] Y (-1)^ 


k=0 


{Na + Nb - Ny{NA - I + ky 
k\{NB - N - ky{NA - ly 


^NB-N-k 


Comparing this result with eq. ()2.8I) derived in [TJ] we can read off the following identity for Rummer’s 
confluent hypergeometric function 


b—a—1 

iFi{a-,b;x)= Y 

j=0 


(-ir(6-l)!(a-l+j)! 1 

j\{b — a — 1 — j)!(a — 1)! x“+-5' 


I X ip' {b - iy{b - a -1+jy 1 

^ - I - jy{b - a - ly x^-^+N 

(A.14) 


for b > a positive integers. This expression seems to be divergent at x = 0. However one can check 
term by term that all negative powers in x cancel in both sums. The existence of such expressions is 
well known for a — b integer or a a positive integer, cf. [32], and can be derived alternatively using the 
recursion for iTi(a;6;x) together with the known initial conditions for iTi(a;a;x) = exp[2;] and for 
iFi(a;a + l;x) in terms of the incomplete Gamma function. 


B Extension of Andreief’s Integration Formula 


For completeness we quote here the generalisation of Andreief’s integration formula derived in [36j . 
as it is used several times in the main text. Let Rj{x), 1 < j < N + k, and Sj{x), 1 < j < N + I, 
be suitable integrable functions and {rab}i<{Yk^^ constant matrices. Then 
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the following integral identity holds, 



= (-l)^'lV!det 


' Rh(x ■ 

det 

Sh(x 

1 l^b<N^k 


n<b<N+l 





J axna\X)ShyX)\^^^^j^_^^ ^^^\l<a<N+k 


|l<b<7V+Z „ 


(B.l) 


In |36] complex integrals were considered while here we restricted the integration to real domains. 
Indeed one can replace the integrals by any linear functionals Rj acting on the functions Sj since 
eq. (IB.lh is only an algebraic identity. It needs the linearity of the integral, the multi-linearity and 
the skew-symmetry of the determinant. Hence the identity (IB.ID would also be true for sums or more 
complicated, in particular higher-dimensional integrals. 


C Expectation Value of Ratios of Characteristic Polynomials 

The approach to calculate Zq^p{X), see eq. ()4.2p . works almost identically as the one in calculating the 
average of the single characteristic polynomial P]y{x), see eq. (I3.15D . The crucial modification enters 
due to the additional characteristic polynomial in the denominator. The source variables yi need an 
imaginary increment to regularise the integral. Thus, let us define L = diag(sign Im(yi),..., sign Im(?/g)) 
that contains the signs of the imaginary parts of the source terms yj. We assume that we have 
positive imaginary increments and g_ negative ones, + q- = q. Then the supersymmetric group 
involved in the partition function Zq\p{X) is U(g+,g_|p), see [53II55] . Indeed we find this group again 
after mapping the average 

Zg|p(X) = J[dA] j[dB]VAiA)VBiB) Sdet-^ (ijv O X - {AA^ + BB^) ® 1^,^) (C.l) 

to superspace. We already rewrote the ratio of characteristic polynomials into the compact notation 
of a superdeterminant to be defined below, as a tensor oi N x N ordinary matrices and of {q\p) x {q\p) 
supermatrices which we will introduce next. 

C.l Brief introduction to the superalgebra of supermatrices 

Let us briefly recall the crucial objects of super analysis and superalgebra with supermatrices and 

introduce our notation. A more thorough introduction into this topic can be found in m- Any 

complex rectangular supermatrix a of superdimensions {qi\pi) x (92IP2) can be arranged into four 
matrix blocks, 

a = f ^bb ctbf \ 2) 

V < 7 FB CTfF / 

The qi x 52 boson-boson block g-qb and the pi x p 2 fermion-fermion block cff comprise commuting 
variables, only. Note that this does not imply that they do not contain nilpotent terms, they indeed 
can. In contrast to the diagonal blocks the qi x p 2 boson-fermion block cjbf and the pi x q 2 fermion- 
boson block (Tfb only consist of anti-commuting and, thus, nilpotent matrix entries. We denote by the 
set Gl(gi|pi; ^21^2) those supermatrices where the matrix entries of ctbb and uff are ordinary complex 
numbers and the matrix entries of cjbf and ctfb are independent complex Grassmann variables. In the 
case that q = qi = q 2 and p = pi = p 2 we abbreviate this set by Gl(g|p). Another important set we 
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need is the coset Herm+(g+, g_|p) = Gl((;+ + g_|p)/U(g+, ( 7 _|p). A supermatrix U G Herm+(g+, g_|p) 
has a boson-boson block which is L-Hermitian, i.e. = LUbbL, and fulfils the positivity condition 
LUbb > 0. We emphasise that these two properties imply that Ubb bas q+ positive real eigenvalues and 
q- negative real eigenvalues. The fermion-fermion block of U is unitary, i.e. Upp = Upp and the off- 
diagonal blocks consist of independent complex Grassmann variables with the condition C/pg = Ubf- 
q- =0 we write Herm_|_(g|p) = Gl(g|p)/U(g|p). A matrix V in the non-compact supergroup 
U(q+,q_|p) satisfies the identity 

V-^ = LV^L, (C.3) 

with the diagonal supermatrix L = diag(L, Ip). 

Two important functions of a supermatrix a G Gl(g|p) are the supertrace, 

StrcT = Tr cjbb — Tr cjff , (C.4) 


and the superdeterminant. 


Sdet((j) = 


det [ctbb — O'BFO'FFf^FB] 
det [uff] 


det [cjbb] 


det [cjff — o'fbo'bb^bf] 


(C.5) 


where here we need cjff and cjbb to be invertible. The definitions are chosen in such a way that 
the properties of cyclic permutation invariance (StrAi? = Sir BA), or factorisation (SdetAi? = 
Sdet ASdet S), and the relation InSdet(A) = STrlnA are natural generalisations from the ordinary 
trace and determinant. Note that for the invariance under cyclic permutation the supermatrices A 
and B can also be rectangular while for the other properties we need square supermatrices. 


C.2 Mapping to superspace 

To keep the calculation simple we omit the normalisation constant of Zq^p{X). In the end we fix it via 
the asymptotic 

lim Z.\JeX) = Sdet-^(A). (G.6) 

In a first step we rewrite the superdeterminant in eq. (|G.ID as a Gaussian integral of a rectangular 
supermatrix V G G1(A^|0; c/|p), i.e. 




jj-^y]e*Tr VLxyt-iTr (AA^+bb^)vlv^ 


/[dyje-Tr^vt 


(C.7) 


The prefactor comes from the factor iL which has to be introduced into the superdeterminant to 
guarantee the positive definiteness of the Hermitian part of the matrix in the Gaussian integral. In 
this way the imaginary increment carries over to superspace. In particular, it allows us to interchange 
the integrals over A and B with those over V since all integrals are absolutely integrable. 

The integrals over A and B are purely Gaussian now, cf. eqs. (12.11) and (IC.7p . such that we can 
perform them and find 


Zgip{X) (X J [dV] exp[iTrVLXV^ def^^ 


+ iVLV^ + iVLV^ 


^B 


(G.8) 


The two determinants are immediately regularised by the imaginary unit in front of VLV^. 
In the next step we apply the duality relation 


det"^^ + iVLV^] = det^^ [Sa] Sdet"^^(lg|p + iLV^J^aV) 


(C.9) 
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and similarly for the determinant with Sb- This relation is true because one can identify det[...] = 
Sdet(...) for boson-boson blocks and use the invariance under cyclic permutation, Sdet(l — YZ) = 
Sdet(l — ZY), reminiscent of the cyclic permutation invariance of the supertrace. Then the partition 
function reads 


Zg\p{X) oc J [dV] exp[iTr yLXy^]Sdef^^(lg|p + iL0s^y)Sdet-^« (l^lp + iLV^J^sV). (C.IO) 
To rewrite the two superdeterminants again as Gaussian integrals, namely as 

/IdlTAlexpl-StrlV'H'^l 


with Wa G G1((7 |p; A^yi|0) and analogously for the one with Sb, we need to discuss if the Hermitian 
part of the boson-boson block of lq\p+iLV^Z,AV is positive definite, in particular if the Hermitian part 
of Ig -|-iLHggS/iHBB is positive definite. Since I^b is an ordinary complex rectangular q x Na matrix 
we need to know if the matrix K = LV^^ZaVbb has complex eigenvalues. Thereby we note that LK 
is Hermitian and positive definite. The L-Hermiticity of K tells us that the eigenvalues of K are either 
real or complex conjugate parts. Moreover, we can block-diagonalise K = U~^'EU by a non-compact 
unitary matrix U G U((/+,g_), i.e. 17~^ = LU^L. In the simplest case of L = diag(lg^, —lg_) the 
matrix H is of the form 



/ diag(xi7,...,xpg^_;) 

0 

0 


0 

diag(x 2 ,i,.. 
-diag(7/2,i,. 

■,X 2 ,i) diag(7/2,i,. 

..,7/2,/) diag(x 2 ,i,. 

• ,y2,i) 
■,X2,l) 

0 


V 0 

0 

diag(x37,.. .,X3^q_-i) 


(C.12) 

with xij,X 2 ,j,X 3 j, 7/2,J G 1^ and 1 = 1,..., [min(( 7 _|_, (?-)/ 2 j . The floor function [. . .J yields the largest 
integer smaller than or equal to its argument. The block diagonalisation with the structure ()G.12D was 
also discussed in m where the situation L = 75 = diag(ln, —In+u) was considered. The situation of a 
general L is related to the structure ()G.12D by a simple permutation of rows and columns such that the 
assumption L = diag(lg_p, —lg_) is not a restriction at all. The positivity of LK = LU~^'EU = WL'BU 
carries over to a positivity condition of LH. This implies two things. First, H and, thus K, has no 
complex conjugated pairs of eigenvalues, i.e. I = 0. Second, the eigenvalues xij are positive and 
the eigenvalues x^^j are negative definite. Hence K has a real spectrum and the real part of each 
eigenvalue of Ig -|- iLV^J^AV = C/“^(lg -|- iE)U is equal to 1. This discussion justifies the Gaussian 
integral (|G.lip and renders the integral absolutely integrable. 

Interchanging the integrals over V with those over Wa G Gl(g|p; A^yi|0) and IFb G Gl((/|p; A^b| 0), 
we integrate over V and find 

ZgipiX) oc j[dWA] j < 8 ) X - < 8 ) WaW^ - Sb ® WnW^). 

(C.13) 

In the final step we apply the superbosonisation formula [26 11391 - 142] for Na,Nb > N > q to replace 
WaW\ and WbW^ by Ua, Ub G Herm+(( 7 |p). This yields 


Zg|p(X) OC j dniUA) j d/r([/B)e"®^’'^^“®*’'^®Sdet^^(C/A)Sdet^s([/B) (C.14) 

xSdet ^(Iat ^ X — Tia ® Ua — Tb G) Ub) ■ 
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(C.15) 


The Haar measure of the coset Herm+(g|p) is explicitly given by [3911421157] 

dn{U) = {27ri)-P SdetP-'i{U)[dU], 

where [dU] is again the flat measure, in particular the product of differentials of all independent 
matrix entries. We emphasise that there is no natural normalisation of the Haar measure on the 
coset Herm+(g|p) when pq > 0; namely the volume of the supergroup U(l|l) vanishes due to Cauchy- 
like integration theorems |58H60j . Hence we choose the normalisation by convenience since p contour 
integrals are involved, see the next subsection. 


C.3 Calculation of the normalisation constant 


Due to the asymptotic (1C.61) the proportionality constant in eq. (IC.14p is the inverse of the following 
integral 

C-i = = I dp{UA)e-^^^^^Sdet^^{UA) j dpi{U^^det^^ {Ub) • (C.16) 

Hence the two integrals in Ua and Ub factorise. 

Let n € N. To calculate the integral 

C-^ = J dp(C/)e-®*’^^SdeC(C) = J , (C.17) 

with U G Herm+(( 7 |p), we first split the supermatrix U into four blocks as in eq. (IC.2h . where Cbb £ 
Herm_|_(g) is Hermitian and positive definite and t/pF S U(p) is unitary. 

Employing the second equality (|C.5p for the superdeterminant we shift f/pp —^ Cpp + t/pBCggt/pg. 
The integrals over the complex Grassmann variables comprised in [/pB become Gaussian such that 

C-i= / [dUBB]e-^^^^^dee-^UBB] [ 

./Herm+(5) •'U(p) f27r*jr^ 

The two remaining integrals are well-known from [341142] . The non-compact integral is related to the 
Laguerre ensemble 

[ [dUBB]e-^^^^^dee-‘^[UBB] = ( TT | TT A2({A4) (C.19) 

Juerm+iq) 

q-1 

= 7r^(n - g-I-/)!. 

1=0 


The prefactor in the first line is the volume of the coset U((;)/[U'^(1) x S(g)] and in the second line we 
get an additional factor from the Selberg integral of Laguerre type [34] . The compact integral is an 
integral over the circular unitary ensemble 


/u(p) (2v^^)P ^ 


n 

\i=i 

p-i 

n 

/=0 




fljj |Ap({e*^“})|2 (C.20) 


i=i 


TT 


{n-q + l)\ 


The contour integrals in the second step are of Selberg-type and were computed in |l2]. The combi¬ 
nation of these results for the constants together with eq. ()C.14p yield the result (|4.7p . 
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D Computation of the Spectral Density in the General Setting 


To compute the spectral density for general empirical matrices and we start from eq. (1321) 
with p = q = 1. We choose the following explicit coordinates for the supermatrices 


Ua/b — 


^ / sa/b ^a/b 


Va/b 


■A’Pa/b 


(D.21) 


with sa,sb G 1K+, 4>a^4>b G and two complex Grassmann variables rjA and riB- The measure 

is given by 

e"‘^^/B(ig^^^d(t>A/Bdr]A/Bd'n*A/B 


dpi^A/B) — 


27r 


(D.22) 


To integrate over the Grassmann variables we have to expand the superdeterminants involved in 
eq. ()4.7I) which is 


Sdet^UA/B = 


sa/b + e 


■A’Pa/i 


^V*a/ b^a/b \^ 


= 1 + 


sa/b^ 


i<t>A/BdAIBVAIB 1 J (D.23) 


and 


Sdet ^(Iat (g) diag(?/,x) - Da ® - Ss (g) C/b) 

= (i^ y Tr (F,S,FA)nln^ 

detfelA, - - sbEb) 1 ^ 

+ ^Tr (F^SATsEA)Tr {F^T.bFs'^b) - Tr {F^'FaFs'^b)^^^ {F^'FbFs'^a) 

-Tr {F^T,aFsT.aF^J:bFsT,b) + Tr 

Here we have introduced the short-hand notation 

F^ = {xIn — e^'^^T.A — and Fg = {ylN — sa'Fa — sb'^bY^- 

We insert these expressions into eq. ()4.7p and integrate over the Grassmann variables, 


(D.24) 


(D.25) 


2tt 2-71 Jq Jq 

det(xljv - F'^^^a - F^^T^b) / - Ti (F^J^aFs^a) + Tr (F^J^bFs^b) 

det(ylAr — saSa — sbSb) \SASBe*‘^^e*'^s sbF^^ ^ saF^^ ^ 


+ (^Tr (T^SaTsSa)^ {F^T,bFsT,b) - Tr (T^SaTsSb)^ (T^SbFsSa) 
—Tr (F^SaTsSaT^SsTjSb) -|- Tr (F^SaT^SsT^SbTjSa)^ 


(D.26) 


The integration over the two angles can be summarised in terms of the expression 



„-iUA(l>A 

2tt 


^—iUBCpB 


exp[e*'^^ + e*"^^] det(xl n - 


-e*‘^«Ss + 5) 


1 / C 

uaIi^b'- V ^ dpAj V ^ dpB 


^B 

det(xlAr - /TaSa - PB^B + S) 


^A/B—0 


(D.27) 
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where S' is a fixed N x N matrix and i'a,i^b S Nq. Then the generating function (jD.26p can be written 
as 


Zi\i{y,x) = 


1 


1 + 


(Na - l)!(iVB - 1)! 

Q \ Na — I / Q —1 


dy-Aj 


1 + 


S/is) 


/•CO rQ. 

/ Asa 
Jo Jo 


:^-Sa-Sb 


dSB 


def^iylN - - sb^b) 


X 


1 


sasb 


1 + 


d 

duA 


1 + 


d 


dfiB 


+ — 1 1 + 
SB 


d \ ^ d‘^ 

d/jB) dzAdz\ SA V duAj OzbOz 


B 




xIn - RaT^A - hB^^B 

V^zaT^a V^zbTjb 

) det 


ylN - saTa - sbT^b ^ ^ i \ 

' dzAdz\dzBdz*^^ 


'\/2z'^Tib 

2 i)J 


(D.28) 


P-A/B—^A/B—^ 


Here we rewrote the integration over the two angles and the Grassmann variables as derivatives with 
respect to the source variables y-A/B za/b- Taking the derivative with respect to x, setting x = y 
and taking the imaginary part in the limit Im(y) —)• 0+ yields the result ()4.10l) . 
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